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Abstract 

Dominant energy subspaces of statistical systems are defined with the help of restrictive condi- 
tions on various characteristics of the energy distribution, such as the probability density and the 
forth order Binder's cumulant. Our analysis generalizes the ideas of the critical minimum energy- 
subspace (CRMES) technique, applied previously to study the specific heat's finite-size scaling. 
Here, we illustrate alternatives that are useful for the analysis of further finite-size anomalies and 
the behavior of the corresponding dominant subspaces is presented for the 2D Baxter- Wu, the 2D 
and 3D Ising models. In order to show that a CRMES technique is adequate for the study of 
magnetic anomalies, we study and test simple methods which provide the means for an accurate 
determination of the energy - order-parameter {E, M) histograms via Wang-Landau random walks. 
The 2D Ising model is used as a test case and it is shown that high-level Wang-Landau sampling 
schemes yield excellent estimates for all magnetic properties. Our estimates compare very well with 
those of the traditional Metropolis method. The relevant dominant energy subspaces and dominant 
magnetization subspaces scale as expected with exponents a/v and j/v respectively. Using the 
Metropolis method we examine the time evolution of the corresponding dominant magnetization 
subspaces and we uncover the reasons behind the inadequacy of the Metropolis method to produce 
a reliable estimation scheme for the tail-regime of the order parameter distribution. 
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I. INTRODUCTION 



Computer simulations based on Monte Carlo sampling methods have increased dramati- 
cally our understanding of the behavior of the standard classical statistical mechanics sys- 
tems (for instance Ising-like models), but also of the more complex systems, such us disor- 
dered media, polymeric and glassy materials. Our main approach, in the past half-century, 
was based on importance sampling in the canonical ensemble. The Metropolis method and 
its variants were, for many years, the main tools in condensed matter physics, particularly 
for the study of critical phenomena However, for complex systems effec- 

tive potentials have a complicated rugged landscape with many minima and maxima which 
become more pronounced with increasing system size. In many such cases the Metropolis 
method and its variants are very inefficient methods. Entropic sampling methods (ESM) are 
alternatives to the importance sampling methods, which at least in principle do not suffer 
from such problems. Of course, for second order phase transitions in unfrustrated systems 
cluster algorithms are quite efficient and have essentially solved the "critical slowing down" 
problem. The performance limitations of flat-histogram or entropic methods have recently 
attracted considerable interest. Even for simple systems, such as the Ising model, such 
methods have tunnelling times in energy space that are not proportional to iV 2 as should be 
expected for a pure random walk but they are proportional to a higher power. Moreover, 
it has been shown that tunnelling times may be strongly depended on the model under 
consideration jf|. Furthermore, in order to apply an ESM the density of states (DOS) of 
the system should be known. 

Over the last decade, several efficient methods that directly calculate the DOS of classical 
statistical models have been developed. A few remarkable examples are the entropic 0, Bl , 
multicanonical £|, broad histogram P, Q], transition matrix h| and Wang-Landau [j2j | 
methods. Using these methods, it is now possible to accurately estimate the DOS of quite 
large classical statistical models Q]. Since for complex systems the Metropolis method may 
be trapped for very long times in non-representative energy subspaces it is reasonable to 
consider an ESM program using an approximate DOS as an alternative to the traditional im- 
portance sampling methods. Although this idea exists in the literature j^J , the effectiveness 
of the various possible implementations have not been exploited by systematic comparative 
studies. For instance, the following two-run strategy may be used in an ESM program. 
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In the first stage an accurate estimation of the DOS of the system under consideration is 
achieved with the help of say the Wang-Landau (WL) method, and in the second stage the 
derived DOS is used in an ESM to estimate further properties of the system, such as the 
order-parameter distribution. In such a two-stage program, the CRMES method recently 



developed by Malakis et al. 



14 I5I ] may be used to restrict the energy space and make such 



a program more efficient for the estimation of the critical behavior of any statistical system. 
Our first objective in this paper is to extent the CRMES condition and to observe to what 
degree appropriately defined energy subspaces are sufficiently large for the estimation of all 
finite-size anomalies. Our second objective focuses on the possibility to obtain all critical 
properties of the system by using a one-run strategy of an ESM, based on a WL random 
walk in a restricted energy subspace. 

The rest of the paper is organized as follows. In Sec. |H] we provide a brief outline of the 
CRMES restriction and give alternative definitions of the dominant energy subspaces. The 
scaling of the extensions of these subspaces is illustrated for the Baxter- Wu and Ising models. 
In Sec. IIII Al we discuss how, by employing the WL and the N-fold implementation of the 
WL scheme in a restricted energy subspace, we may generate approximations of the DOS of 
the system and at the same time obtain energy - order-parameter (E, M) histograms. It is 
suggested that this one-run WL strategy yields good estimates of all magnetic properties. 
In Sec. IIII Bl we consider the 2D Ising model as a test case. Our estimates are compared 
with those of the traditional Metropolis method and the expected magnetic scaling behavior 
is recovered. Finally, in Sec. IIII CI we study our dominant subspace method in the energy 
and the order-parameter space. The subspaces sufficient for an accurate estimation of mag- 
netic finite-size anomalies are determined and their scaling is analyzed. The tail-regime of 
the order-parameter critical distribution is briefly discussed and the shortcomings of the 
Metropolis method are clarified. Our conclusions are summarized in Sec. IIVI 



II. ALTERNATIVES FOR RESTRICTING THE ENERGY SPACE 



Let us start by recalling the original CRMES restriction of Malakis et al. [14 1. Assume 
that E denotes the value of energy producing the maximum term in the partition function 
of the statistical model, for instance the Ising model, at some temperature of interest. Since 
we deal with a finite system of linear size L, we are interested in the properties (finite-size 
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anomalies) near some pseudocritical temperature T£, which in general depends on L but also 
on the property studied. Thus, for the specific heat peaks let us use the notation T£[C] = 
T L)C for the pseudocritical temperature and define a set of approximations by restricting the 
statistical sums to energy subranges around the value E = i£(T£[C]) = E(Tl^c)- Let these 
subranges of the total energy range (E min , E max ) be denoted as follows: 

(E-,E+), E± = E±A ± , A ± >0 (la) 

Accordingly, the peaks are approximated by: 

C L {E_,E + ) = C L (A ± ) = N- l T- 2 <{ Z" 1 J] £ 2 exp [$(£)] - Z' 1 ^ Eexp [$(£)] 

(lb) 

E+ 

${E) = \S{E) - (3E\ - [S(E) -0e\, Z = Y, e*P ME)] (lc) 

e_ 

Since by definition <&(E) is negative we can easily see that for large lattices extreme values 
of energy (far from E) will have an extremely small contribution to the statistical sums, 
since these terms decrease exponentially fast with the distance from E. It follows that, if 
we request a specified accuracy, then we may restrict the necessary energy range in which 
DOS should be sampled. To introduce the minimum energy subspace (MES), we impose 
the condition: 

|C£(A ± ) 



CI 



- 1 



< r (2) 



where r measures the relative error and it will be set equal to a small number (r = 10 -6 ), and 
C* L is the value of the maximum of the specific heat obtained by using the total energy range. 
With the help of a convenient definition, we can specify the minimum energy subspaces 
satisfying the above condition. Their finite-size extensions will be denoted by: 

(AE) C = (AE) Chr ee mm{E + - EJ) (3) 

Demanding the same level of accuracy (r) for all lattice sizes, we produce a size-dependence 
on all parameters of the above energy ranges and in particular the extensions of the critical 
MES should obey the scaling law [14]: 

(AE)r 

*C EE i-^ « L~ (4) 



the MES we may follow successive minimal approxi- 



In order to determine the location o 
mations to the specific heat peak [14]: 



C L (j) = C L (Aj,A+), Af +1 = Af ±Sf +1 , A± = 0, j = l,2,... 



(5a) 



where one the above ^-increments is chosen to be 1 and the other according to which side 
of E is producing at the current stage the best approximation: 

(0+ +1 = l, ej +1 = 0) <*\C L - C L (AJ, At + 1) [<| - C L (AJ + 1, At) | 

= 0, ej +1 = 1) <*\C L - C L (Aj, A+ + 1) |>| C L - C L {Aj + 1, A+) | (5b) 

The above defines a sequence of relative errors for the specific heat peaks (rj): 

ClU) 



c* 



1 



(5c) 



and the MES is the subspace centered at E corresponding to the first member of the above 
sequence (J3J) satisfying: r 3 - < r. The location of these subspaces can be predicted either by 
extrapolation, from smaller lattices, or by using the early-stage DOS approximation of the 
WL method |lj|. Using a sufficiently wider subspace we can also accurately estimate their 
extensions and verify the scaling law (0J) as shown in Ref. [l^ . 



The above scheme has been tested for the Ising and Baxter- Wu models [1J, [15j and it 
has been shown that this particular rule provides very accurate extensions satisfying quite 
well the expected scaling law (j3J). However, one may conceive other ways for specifying the 
locations of the energy subspaces that will essentially produce comparable approximations. 
A simple idea is to use a condition based on the energy probability density (fTT,c{E) oc $(E)) 
meaning the application of Eq. at a particular pseudocritical temperature T^c- That is, 
we may define the end-points (E±) of the subspaces by simply comparing the corresponding 
probability densities with the maximum at the energy E: 



E±: exp{$(£ ± )}<r 



(6) 



Applying this restriction at a particular (pseudocritical) temperature (T[) is much easier 
than applying the successive minimal approximations described in Eq. © and it will pro- 
duce comparable approximations for the specific heat maxima. The corresponding scaling 
variable ^f T * c <e) for the resulting MES should also obey the law (J3J). Another interesting 
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pseudocritical temperature, is the temperature corresponding to the forth order Binder's 
cumulant for the energy distribution. This reduced cumulant is defined by: 

(£ 4 > 



V L (T) = 1 



3(E 2 y 



(7) 



and it is well known jisl . 16, ^| that this quantity has a minimum, V£, at a pseudocritical 
temperature T£[V] = T^y. In the thermodynamic limit this temperature also approaches 
the critical temperature but for finite lattices is different from the pseudocritical temperature 
Tl,c of the specific heat. Thus, the maximum term of the partition function, corresponding 
to the temperature of the cumulant finite-size anomaly, will be now located at a different 
value of the energy spectrum (say at E(Tl v))- Therefore, if we follow the probability 
density criterion described above in Eq. (JHJ) to define the CRMES at these temperatures we 
will generally find subspaces that do not coincide with the subspaces for the specific heat. 
However, for large lattices the non-overlapping parts of these subspaces are relatively very 
small and we may run the WL algorithm in the union of these subspaces in order to study 
both properties. 

Let us now, follow an analogous approach with the successive minimal approximations of 
Eq. (J3J) for the cumulant anomaly. It is of interest to note that if we define the CRMES by 
a similar condition: 



V? 



< r 



(8) 



then, for sufficiently large lattices, the scaling law (0J) for the corresponding extensions 
(we may now use the = L J V variable) will not be obeyed. This is due to the fact 
that the cumulant goes asymptotically to a finite value (V^ = V^, = 2/3) and therefore the 
relative accuracy criterion (JHJ) is not appropriate for large lattices or, in other words, becomes 
ineffective (see also the discussion bellow). We may now introduce the finite-size distance 
from the asymptotic value of the Binder's parameter, (V^> — V£), to our considerations and 
modify the criterion (JHJ) as follows: 

VZ(A±) 



YL 
v* 



< r 



(9) 



It appears (see bellow) that this option makes the extensions of the resulting subspaces 
to follow very well the scaling law (0J) and for the 2D Ising model we obtain an almost 
perfect coincidence with the extensions of the corresponding subspaces obtained for the 



specific heat. The corresponding scaling variable will be denoted by ^ 



(AJ5)? 



For 



a second-order transition the limiting value of the energy cumulant is known (V^ = 2/3) 
and this makes the implementation of the scheme © possible. 

The above alternative definitions for the CRMES were applied to the 2D Ising model 
using the DOS data of Ref. jlj (L = 10 - 100) but also new data up to L = 200 (L = 
120, 140, 160, 180,200). Fig. Q illustrates the scaling behavior of the corresponding scaling 
variables for the CRMES of the specific heat at its pseudocritical temperature T^c using the 
minimum subspaces resulting from Eq. (J5cj) . and also the minimum subspaces resulting, at 
this temperature, from the probability density condition studied in Eq. (JBJ). The same figure 
displays the estimates for the two options (jHJ) and (JHJ) for the CRMES corresponding to the 
minimum of the Binder's parameter at its pseudocritical temperature T^y. The expected 
logarithmic scaling law [lJ . Irsj ] is obeyed well for all definitions with a clear exception of 
the cumulant condition (jHJ) in which case the scaling variable shows a clear decline from 
the logarithmic divergence for large lattices. This is due to the fact that this parameter 
approaches the finite value = 2/3, and the condition defined in Eq. (jHJ) becomes ineffective 
for large lattices. Note that with increasing lattice size several significant figures of this 
finite number are determined from a small part of the dominant energy subspace. The trick 
proposed in Eq. Q not only keeps the scheme in the proper scaling form but also it is 
remarkable that the resulting extensions of the CRMES of the specific heat and the Binder's 
parameter completely coincide for the 2D Ising model. This coincidence occurs also for the 
Baxter- Wu model, but not for the 3D Ising model (see bellow). Note however, that the 
corresponding pseudocritical temperatures are different and their CRMES do not coincide 
(for L = 100 their displacement is 20 energy levels), only their extensions are equal. Fig. [21 
illustrates the behavior of the same scaling variables for the 3D Ising model using the DOS 
data of Ref. The situation is very similar and is described by the power law (J3J) as 

discussed in Ref. jl^j]. Noteworthy, that the cumulant condition (jHJ) leads to a clear decline 
from the appropriate power law for large lattices and the trick proposed in Eq. (JHJ) seems 
to yield the expected critical behavior. Finall y F ig. |H1 presents the analogous graphs for 
the Baxter- Wu model using the DOS data of Again the scaling variables appear to 
follow the expected scaling law (JU) and a decline in the case of the condition (JHJ) is observed. 
However, this decline is weaker for the Baxter- Wu model due to the fact that the cumulant 
minimum is quite deeper for this model, that is the difference (V^ — V£) is relatively larger 
for the Baxter- Wu model. Fitting these data to the expected power law (see the discussion 
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in the caption of Fig. El) we find that the best estimate for the exponent ajv is produced 
from the original [l4| , Il5j | minimum subspaces of the specific heat. As mentioned above 
the extensions of the CRMES determined for the cumulant by the condition Q for the 2D 
Baxter- Wu model coincide with the extensions of the CRMES of the specific heat (Eq. (0)) 
of the model. 

III. ENTROPIC SAMPLING VIA WL RANDOM WALKS 
A. Microcanonical estimators via the WL scheme 

As mentioned in the introduction, the ESM method using an exact (if available) or an 
accurate approximation of the DOS may be considered as an alternative to the various 
importance sampling methods used in the literature to estimate canonical averages and 
probability distributions of macroscopic thermodynamic variables. Here we shall examine 
the idea of producing accurate estimates for finite-size magnetic anomalies by using a simple 
ESM based on the WL random walk in an appropriately restricted energy subspace (Ei, E 2 ). 
We shall also test our results by comparing to the Metropolis method l|. 

We implement a WL random walk in a restricted energy subspace (Ei,E 2 ) and at the 
same time we accumulate data for the two-parameter (E, M) histogram. For a large system 



we employ a multi-range algorithm |12j in which independent random walks are used for 
different energy subintervals and the resultant pieces are then combined to obtain the DOS 
in (Ei,E 2 ). The WL modification factor (fj) is reduced at the jth iteration according to: 
fi = e, f) — > 7 = 2 J f in and for the histogram flatness we use a 5% criterion as in 



our previous studies [1J, [la, |19( . Note that, the detailed balance condition is satisfied in the 
limit / — > 1. Let the exact density of states be denoted by G(E) and the DOS of the above 
WL process be denoted by Gwl{E). Similarly, let Hwl(E, M) be the histogram produced 
during the WL process by a specific recipe which will be described below. 

The resulting approximation of the DOS and the corresponding E, M histograms may 
be used to estimate the magnetic properties of the system in a temperature range, which 
is covered, by the restricted energy subspace (Ei,E 2 ). Canonical averages will be then 
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approximated as follows: 

(M n, _ E E (M n ) E G(E)e-^ EE, iEu E 2) (M n )E,WLG WL (E)e-^ 

{ } E E G(E)e^ E mEl ,E 2) Gw L (E)e-^ { > 

where the micro canonical averages {M n ) E are obtained from the Hwl{E, M) histograms as: 

(M n ) E S (M n ) E>WL = V M n HwL{E ^\ H WL {E) = V H WL (E, M) (10b) 

and the summation in M runs over all values generated during the process in the restricted 
energy subspace (Ei,E 2 ). 

The accuracy of the magnetic properties obtained from the above averaging process will 
depend on many factors. Firstly, the used energy subspace restricts the temperature range 
for which such approximations may be accurate. This restriction has as a result that the 
process will not visit all possible values of M, but this fact is of no consequence for the 
accuracy of the magnetic properties at the temperature range of interest, as far as the es- 
timated DOS is accurate. Secondly, the accuracy of the above microcanonical estimators 
will, as usually, depend on the total number of visits to a given energy level (H WL (E)), and 
also to the number of different spin states visited within this energy level. However, these 
are statistical fluctuations inherent in any Monte Carlo method and we should expect im- 
provement by increasing the number of repetitions of the process. At this point let us point 
out that statistical fluctuations may be reduced, as usually, by multiple measurements but 



also by using some refinements of the original WL algorithm [20, |2l| . An illustrative figure 
including such a refinement will be presented in the next subsection. Finally, the construc- 
tion of reliable (uniform) approximations for microcanonical averages is an important open 
problem, discussed recently in some detail by P. M. C. de Oliveira |22J]. The microcanonical 
approximations ((M n ) EtWL ), appearing in Eq. (JTUJ), and used in this paper are obtained 
from the WL multi-range process, and the corresponding N-fold version of Shulz et al. j3| . 
The recipes employed are outlined and tested in the sequel. 

There are mainly three categories of microcanonical simulation approaches. In the first, 
one tries to satisfy completely the fixed-energy constraint, a typical example is the Q2R 
cellular automaton [3, l^ . In the second, one tries to mildly relax the energy constraint 
relatively small energy-windows in order to avoid ergodic problems, as done by 
2q in his "demon" method. Finally, the fixed-energy constraint is completely re- 



by usin. 
Creutz 
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laxed 27[ and various ideas have been tried, ranging from the tuned temperature canonical 
approach of Oliveira 



to the flat histogram recipe of J.-S. Wang [2£ 
In order to present some justification for the scheme of Eqs. (|10a|) - ([10b[) let us suppose 
that we know the exact DOS and we are about to perform an infinitely long entropic sampling 
in the restricted energy subspace (Ei,E 2 ). In this case we can write for the Monte Carlo 
process with sampling probability Pi oc 1/G{E) 



mo 

4]: 



(Mn) = E iEl , E2) (M n ) E G(E)e-^ YJLm-P- 1 ^ E {El , E2) (M-)E,ESMG(E)e-^ 
[ } E [El , E2) G(E)e-^ Tl.Pr 1 ^ E (EuE2) G(E)e^ 

(11a) 



The last approximation in Eq. ()lla|) assumes that in the limit of an infinite Markovian 
chain the Hesm(E) histogram is perfectly flat. Accordingly, Hesm(E) has been set equal 
to a constant in the denominator (in replacing the sum over the Af sampled microstates by 
a sum over energies) and then moved to the numerator in order to form the ESM micro- 
canonical averages defined by: 

(M u ) e , E sm = J2 M" ^ A/( ^ M) , H E sm(E) = V H ES m(E, M) (lib) 

The above observation shows that the ESM microcanonical average should be an unbiased 
estimator for the fundamental microcanonical average: 

(M n ) E , ES M^ (M n ) E (12) 

Therefore, Eq. ()12|) provides the essential theoretical support for using a multi-range 
WL process (at its late stages) to obtain microcanonical simulators. It is reasonable to 
expect that the high-level stages (j ^> 1) of the WL process will resemble the dynamics of 
the ESM and therefore will produce good approximations for the microcanonical aver ages . 
This approach is similar in many respects to the flat histogram method of J.-S. Wang [28]. 
Since the resemblance of the WL process with the ESM depends on the value of the control 
parameter (fj), we shall classify our recipes by using the j-range utilized for updating the 
(E, M) histogram during the WL process. Thus, if all accepted microstates of the WL 
process during the j-range (j = Ji n u, J fin) are used, for updating the (E, M) histogram, 
the resulting recipes will be denoted by WL(J init — J/, n ). When using the N-fold version of 
the WL process we always start with several (Jwl) WL j-iterations and then continue the 
process from the next level (JN-foid = J\vl + 1) by carrying out further N-fold WL iterations 
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(j — JN-foid, J fin)- 111 this recipe, we shall use only the N-fold iterations for updating 
the (E, M) histogram and time-weight the hi stog ram by using the life-time of microstates 
calculated according to the N-fold method 0, [ljQ] . The resulting recipes are denoted by 
WL(N-fold)(J*_ /oW - J fin ). 

The above described schemes were tested by using Eq. (II (J j) to obtain certain magnetic 
properties of the 2D Ising model (for lattices with sizes L = 20 — 100) and compare to the 
estimates obtained by the Metropolis method. The finite-size anomaly of the susceptibility 
and its value at the exact critical temperature, as well as the value of the order parameter 
also at the exact critical temperature were used in these tests and appear in the following 
subsection. For the magnetic susceptibility we have used the definition: 

XL(T) = ±(m^Ml) (13) 



N \ T 
and for the order parameter: 

m T =(\M\/N) T (14) 



B. Metropolis versus WL - CRMES schemes. A comparative study 

The estimates of the Metropolis method were obtained as follows. First an initial equi- 
libration period of 50 x L 2 usual Monte Carlo steps (lattice sweeps) was applied without 
updating the histograms. After thermalization, the updating of the histograms was applied 
in every Monte Carlo step, while the magnetic properties, the order-parameter distributions 
and the time evolving dominant M-subspaces were determined and observed in time steps 
of 50 x L 2 Monte Carlo sweeps. The time evolution lasted a total of 300 such time steps 
for all lattice sizes (for L = 70, 120 and L = 140, see also the discussion bellow). Thus, for 
a lattice of linear size L = 100, the above time evolution accounts to a total of 1.5 x 10 8 
lattice sweeps. 

The estimates of the WL multi-range process were obtained using, for each lattice size, 
30 independent random walks in the appropriate (Ex, E 2 ) energy subspace. These subspaces 
were chosen carefully to cover a wide temperature range close to the critical point, so that the 
DOS and the H(E, M) histograms produced would yield accurate estimates of all thermal 
and magnetic properties in the temperature-range. The energy subspaces were wide enough 
to produce relative accuracies , within the scheme, far beyond the criterion r = 10 -6 , which 
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was finally applied to determine the dominant subspaces. The time requirements of the 
described Metropolis estimation were notably greater than the described WL multi-range 
process of 30 independent random walks. For L = 100 the WL scheme was about 3 times 
faster than the Metropolis scheme. Note that the Metropolis method produces estimates 
only for one particular temperature and not for a wide temperature-range. For a lattice 
of linear size L = 140 the above described WL scheme is at least 10 times faster than the 
described Metropolis scheme. 

The estimates derived from the WL recipes appear to have relatively small deviations 
from the corresponding Metropolis estimates and, with a notable exception, they seem to 
be within the given error bounds. The comparison with the Metropolis method is presented 
in Figs. 14151 and lol The relative variations shown in these graphs are defined with respect to 
the Metropolis estimates by assuming that these are more accurate, i.e we define: 

= n ( 15 ) 

VMeir 

where Q = itlt and/or Q = Xt- The error bounds used in these graphs refer to the 
Metropolis estimates in the observed equilibrium regime which is roughly defined as the 
last almost flat part (t = 150 — 300) in the above described 300 time-steps. To define the 
error bounds we have used a confidence level of 5 standard deviations obtained in this wide 
time- window. 

Fig. |3] shows the relative deviations from the Metropolis mean values of the order- 
parameter at the exact critical temperature for four recipes of the WL scheme, as indicated 
in this figure. The Metropolis estimates have very small error bounds of the order of 0.5% or 
less, besides the fact that we have applied the above mentioned demanding confidence limit. 
The deviations of the WL schemes are reasonable (of the order of 1% or less) with the clear 
exception of the case WL(l-24) in which the histogram's (H(E,M)) updating started from 
the very early stage of the WL process. For large lattices, this recipe seems to produce a 
significant overestimation of the order-parameter, at the exact critical temperature, and this 
is enhanced with increasing lattice size. Since, the detailed balance is strongly violated at 
the early stages of the WL process, the observed failure of this recipe should be attributed 
to the detailed balance violation. The related overestimation may be possibly a result of an 
oversampling distortion of large magnetization values at the low energy part of the energy 
range used. Distortions stemming from the violation of the detailed balance condition are 
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difficult to explain and in general their origin is a subtle matter. However, such systematic 
distortions are not observed for the other three recipes appearing in Fig. 0] We consider this 
as a first strong indication that the weak violation of the detailed balance condition for these 
high-level WL schemes is not statistically significant. The behavior of the relative deviations 
for the susceptibilities at the exact critical temperature in Fig. and at the pseudocritical 
temperatures in Fig. |U] appear in general, much better from those shown in Fig. |3J Again 
the distortions of the WL(l-24) scheme become pronounced as we increase the lattice size. 

At this point let us try to observe in more detail the effect of the WL iteration level 
on an estimated magnetic property and also the effect of one of the simplest refinements 
of the WL algorithm for the square Ising lattice of size L = 30. Using an accurate DOS, 
Gwl, obtained from a previous WL run = 24), we have calculated, with the help 

of Eq. ([lip, the critical susceptibility xt c as a function of the WL iteration level in a new 
WL diffusion process. In this new multi-range (mr) process each WL iteration level is 
repeated 30 times for each energy subinterval. Fig. [7| shows the evolution of the critical 
susceptibility to its equilibrium value for five cases. In the first case (mr-WLl) the updating 
of (E, M) histograms follows the recipe WL(1 — 24) and in the second case (mr-WL2) the 
recipe WL(12 — 24). The third case (mr-WL2iS) follows a refinement of the WL algorithm 
proposed by Zhou and Bhatt [20]. The additional element of the algorithm is that now we 
allow for a number S (S = 16) spin-flips between successive records in the histograms (and in 
the DOS modification). Introduction of such a separation diminishes systematic errors due 
to the correlation between adjacent records as shown in Ref. j^J. From the first case we note 
that starting the (E, M) histogram updating process at the early stage of the WL difusion is 
analogous to adding a "non-equilibrium memory effect" in our cumulative histograms. This 
early effect is stronger when the WL algorithm is used in a simple one-range fashion, as our 
simulations have shown. It is also apparent that the refinement introduced by the separation 
S clearly improves, in the cost of the additional spin-flips, the behavior of the evolution of 
the magnetic susceptibility \t c towards its final equilibrium value, which otherwise (5 = 1) 
is obtained in a longer run. 

Zhou and Bhatt [20] have given a proof of the convergence of the WL algorithm and found 
that the fluctuation of the histogram is proportional to 1/ '\/IHT where / is the modification 
factor. This has been recently confirmed by numerical tests |2l| . where it was shown that 
the criterion for reducing / goes beyond the "flat histogram" idea and that "enough statis- 
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tics" should be obtained in each WL iteration. According to Lee, Okabe, and Landau [21 1 
an optimal algorithm will stop the simulation as soon as the histogram fluctuation at the 
jth iteration, denoted as AHj = J2 E (Hj(E) — min E {Hj(E)}), becomes saturated. In order 
to observe whether we will have a noticeable improvement on the magnetic susceptibility 
estimates by applying the proposed entropic sampling scheme (WL( 12-24)) for longer sim- 
ulation times, at each j{= 12 — 24) iteration, we present also in Fig. [7| the cases WL2S(tj) 
and WL2iS(3tj). Both these runs were performed via a simple one-range WL scheme, using 
at each iteration more Monte Carlo steps per spin (MCSS) than required for the saturation 
of AHj. Moreover, the latter case uses three times the number of MCSS of the former. 
The inset of Fig. [7| illustrates the clear saturation of AHj =12 for these runs (the solid (dot- 
ted) line corresponds to duration tj {3tj)). Comparing the estimates of the last three cases 
shown in Fig. [7| (S = 16) we may draw the following conclusions. Firstly, the multi-range 
approach applied in the case mr-WL2S gives comparable estimates to those obtained from 
the one-range implementation of the WL scheme and secondly, a much longer run, such as 
the WL2iS(3t,), does not yield a noticeable improvement. It appears that an optimal and 
quite accurate implementation of the proposed entropic scheme can be constructed using 
a multi-range H(E, M) histogram updating process, during the high-level "well saturated" 
WL iterations. 

The recent combination of Lee's entropic sampling with the WL algorithm presented in 



Ref. 21( may be also implemented to test and possibly improve the proposed here CRMES 
entropic scheme. Finally, the adaptive algorithm of Trebst et al. [29] is of particular interest 
for a CRMES implementation and would be also stimulating to compare the "bottleneck re- 
gion" , or region of minimum diffusivity, of this method with the dominant energy subspaces 
as defined in this manuscript. Therefore, we conclude that, the high-level WL (CRMES) 
schemes are reasonable alternatives to the Metropolis method. The estimates for thermo- 
dynamic parameters involving higher moments of critical distributions appear to be of the 
same or even better accuracy with those corresponding to the traditional method. 

C. Energy and order-parameter dominant subspaces 

The energy E producing the maximum term in the partition function at the pseudocritical 
temperature of the susceptibility T£ [x] = Ti tX may be easily located. Thus, we may apply a 
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minimal approximation scheme analogous to that of Eq. (J3J) to observe the evolution of the 
susceptibility maximum as we expand the energy subspaces centered around E. Now the 
value of the susceptibility is used in place of the specific heat and the resulting CRMES is 
the subspace centered at E corresponding to the first member of the sequence satisfying: 

** (A±) - 1 <r (16) 

Provided that our initial guess for (Ei, E 2 ) is wide enough we also obtain accurate estimates 
for the finite-size extension of these subspaces and we would expect that the relevant scaling 
variable ^f x = Li x would obey the scaling law (J3|. Therefore, we may restate the scaling 
law as: _ 

* eS iM)|«L? (17) 

where 0(= Ql) denotes the finite-size value of some thermodynamic variable close to a crit- 
ical point and (AE)q is the extension for appropriately defined minimum energy subspaces. 
In analogy with our findings jl^ . for a diverging specific heat behavior, we expect that a 
diverging susceptibility and the criterion (|16| will yield extensions scaling according to the 
above law. The alternative method described in Sec. |H] which employs the energy density 
function (see Eq. EJ) may be also used and is easier to implement. 

Fig. |H1 illustrates the scaling behavior of the extension of the CRMES defined with the 
help of the specific heat maxima (Eq. (JHJ) and the corresponding CRMES defined with the 
help of the magnetic susceptibility maxima, as discussed above (Eq. llfijh The corresponding 
scaling variables should be expected to obey the scaling law (fTTj) . and for the 2D Ising model 



the well known logarithmic law 



14. 



la ]. As seen from this figure, this logarithmic law is 



well obeyed for both restrictive schemes as expected. 

Finally, we may describe a procedure for specifying the dominant subspace in the order- 
parameter space. We assume that the energy subspace {E^E^) is sufficiently broad to 
approximate to the desired degree of accuracy the probability density of the order-parameter 
at some temperature of interest by: 

H WL (E,M) n ( -p\ -PE 
p ( M \ ~ ^mEM H WL (E) ^WL{U)e 

Pt{M) = Z mEl , E2) G WL (E)e-^ (18) 

Next, we find the value M that maximizes the above density, at the pseudocritical temper- 
ature T L)X (or some other temperature, for instance the exact critical temperature), and we 
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TABLE I: 3 = ~ aL w , exact w = I = 1.75. Data fitted: L = 50 - 100, T = T c . 





Metropolis 


WL(l-24) 


WL(12-24) 


WL(N-fold: 14-26) 


a 


0.59(6) 


0.31(9) 


0.54(1) 


0.54(1) 


w 


1.73(2) 


1.89(6) 


1.76(1) 


1.76(1) 



locate the end-points (M±) of the magnetic critical subspaces by: 

~ P Tr (M±) , , 

M± : L ' xV ^ < r 19 

~ Pt l JM) 

The above condition is in full analogy with that of Eq. © applied there in the energy 
space and since we will now consider only dominant M-subspaces, defined with the help of 
the above probability density criterion, we shall avoid in our notation the explicit reference 
to the probability density. Accordingly, we denote the extension of the resulting magnetic 
subspaces by: 

(AM) Tlx = (AM) Ptl JM) ee min (M + - M_) (20) 
and we should now examine for a scaling law of the form: 

EE E Ptl JM) EE « L l (21) 

Fig. El illustrates the scaling behavior of the critical minimum magnetic subspaces (Cr- 
MMS) obtained using the magnetic space restriction (|19p. defined above, and the accuracy 
level r = 10 -6 at the pseudocritical temperatures of the susceptibilities. The behavior of 
the high-level WL recipes is very good and provides quite accurate estimates for the critical 
exponent "f/v. The line shown is the power law 0.525 • L 175 which is obtained by fixing the 
exponent to 1.75 and fitting the WL(N- fold: 12-24) data. This is almost identical with the 
power law obtained by using the exponent as a free parameter, which yields 0.535 • L L746 . 
Note that the deviations of the WL(l-24) scheme are not observable for small lattices. How- 
ever, with increasing lattice size (L > 60) they become quite apparent. Fig. [TU] presents 
a similar illustration at the exact critical temperature comparing now not only the WL 
recipes but also the Metropolis method. The deviations from the expected scaling law are 
now quite apparent not only for the "bad" WL(l-24) recipe but also for the Metropolis 
method. Table |l] gives estimates of the critical exponent obtained from the schemes shown 
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in Fig. El using only the intermediate data L = 50 — 100 in which these deviations are still 
moderate. From this table the overestimation of the WL(l-24) scheme but also the tendency 
of underestimation of the Metropolis scheme is quite obvious. The Metropolis data used for 
L = 20 — 100 were obtained from the CrMMS developed at the end of the 300 time-steps 
described in Sec. 1111 Bl It is of some interest to list here some details of the end-points of the 
dominant magnetization subspaces. Let us consider the case L = 100 as an example. The 
broad energy range used in our runs for the WL process was (ie = 850, ie = 2150), where 
the counting of energy states is given by ie = (E + 2iV)/4 + 1. The dominant energy ranges 
are (a) at the pseudocritical temperature of the specific heat (ie = 1061(3), ie = 1947(3)) 
and (b) at the pseudocritical temperature of the susceptibility (ie = 1142(3), ie = 1959(3)). 
The maximum value of the magnetization sampled in the process was \M\/N = 0.912 and 
the minimum value \M\/N = 0. Defining the dominant magnetization space with the help 
of Eq. (|T9"j) at the pseudocritical temperature of the susceptibility, we locate this subspace 
as (\M\/N = 0, |M|/iV = 0.816). This shows that the so defined dominant magnetization 
subspace (CrMMS) is a subset of the sampled values. Note that for all lattice sizes the 
left-end of the dominant magnetization subspace is \M\/N = 0. 

Let us now present the Metropolis time-evolution of the dominant M-subspaces since this 
development offers a didactic example of the very slow tail-convergence of the traditional 
importance sampling methods. Fig. ^2 illustrates the slow equilibration process of the 
Metropolis algorithm at the tails of the order-parameter distribution. To draw this graph 
the time-developing CrMMS, for the two r-levels (r = 10~ 4 and r = 10~ 6 ) shown, was 
divided by the corresponding CrMMS predicted by the high-level WL(N- fold: 12-24) scheme 
which appear to have very small errors for the lattice sizes shown. Thus, considering these 
later CrMMS as exact the Metropolis relative dominant M-subspaces should grow in time 
towards the value 1. For the value r = 10~ 4 this is almost achieved at the time t = 150 for 
both lattices shown, L = 70 and L = 100. However, for r = 10~ 6 , we observe a very slow 
relaxation process which persists even if we increase the observation time. This is obvious 
in Fig. IT2l where, for the smaller lattice L = 70, the time duration of the process has been 
increased up to t = 900. Note that, if one was observing the equilibration process of the 
algorithm with respect to the mean value of the order-parameter, he would then have been 
convinced that equilibrium of this quantity has been attained well before the time t = 150. 

This situation is due to the very slow equilibration at the tails of the distribution and in 
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particular at the large-M tail. In fact the time-expansion of the CrMMS in the Metropolis 
process is due to the gradual movement of the right-end of the magnetization range to larger 
values, since the left-end is \M\/N = from the very first stages of the process. Returning to 
the large-M range is a rare event for the Metropolis algorithm and this makes this traditional 
scheme inaccurate in the far tail-regime, but also very inefficient for the study of the tails 
of the critical distributions. The data points L = 120 and L = 140 for the Metropolis 
algorithm included in Fig. El were obtained by using runs with approximately equal time 
requirements (about 40 hours in a Pentium IV 3GHz) with the WL(N-fold: 14-26) scheme. 
For the case L = 140 this corresponds to only 30 time-steps with the developing relative 
CrMMS (r = 10~ 6 ) having hardy approached the value of 0.97 only. Even by using a much 
longer run (for instance 100 times longer) the Metropolis algorithm will not give an adequate 
description of the far-tail regime. 



IV. CONCLUDING REMARKS 



Critical dominant energy and order-parameter subspaces can be defined by various alter- 
native restrictive schemes, as shown in this work. In this way it is possible to optimize the 
Monte Carlo schemes and study simultaneously all finite-size anomalies of statistical models. 
The finite-size extensions of the dominant energy and order-parameter subspaces scale with 
exponents ajv and 7/1/ respectively. Our experience with this subject leads us to conclude 
that the extensions of these dominant subspaces are more accurate than the estimates of the 
corresponding thermodynamic variables (specific heat and susceptibility), establishing the 
critical minimum subspace method as a new alternative for the estimation of the associated 
critical exponents from finite-size Monte Carlo data. 

The presented clarification and generalization of the CRMES method greatly speeds up 
the Monte Carlo simulations in many applications of the methods determining the spectral 
densities in classical statistical models. Furthermore, the presented one-run high-level WL 
entropic sampling schemes provide efficient alternatives when carried out in appropriately 
defined dominant subspaces. We expect that for complex systems with long trapping times, 
these schemes will appear to be much more advantageous in almost all respects. This last 
expectation has been verified by our studies of the random-field Ising model, which will be 
published shortly. Moreover, these methods have general advantages, the most important 
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of these being: (a) the fact that one can improve the H(E, M) histograms by repeated 
application of the method (at the same time we improve the accuracy of the DOS in the 
energy space), and (b) the fact that their implementation in a sufficiently broad energy 
subspace provides data for calculating all finite-size properties of the statistical system, which 
are relevant for the prediction of the asymptotic critical behavior. Overall, we envisage that 
our study can be further utilized in many ways for the investigation of the critical behavior 
of statistical models in future researches. 
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FIG. 1: 2D Ising model (L = 10 — 100 [l4|, L = 120, 140, ...,200 new data). Behavior of scaling 
variables corresponding to various alternative definitions of dominant energy subspaces defined 
in the text. VP/ T » c [e) is t ne scaling variable defined with the help of the restriction © at the 
pseudocritical temperature Tl,c- Correspondingly, ^>c is the scaling variable defined from Eq. (J5J) 
at the pseudocritical temperature T^c an d ^v-y^ is t ne scaling variable defined from Eq. © at 
the pseudocritical temperature T^y. Finally, tyy is the scaling variable defined from Eq. (jHJ) at the 
pseudocritical temperature T^ y. Note the strong decline of this last variable from the expected 
asymptotic law, as explained in the text. The log scale in x-axis is used in order to facilitate the 
observation of the logarithmic behavior. 
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FIG. 2: 3D Ising model (L = 4 — 32 [141] ). Behavior of the scaling variables \E'c', ^v-,Voo and ^y, 
as in Fig. ^ Note again the clear decline of the variable defined by Eq. (|SJ). The log- log scale 
is used in order to observe the expected power law. 
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FIG. 3: 2D Baxter- Wu model (L = 42,51,60,69,75,78,84,87 and L = 96 Behavior of the 

scaling variables ^/ T * c (E)i ^C, ^v-y^ an d VPy, again as in Fig. ^ The modest decline of the 
variable Vl/y defined in Eq. (JSJ) is due to the fact that the cumulant minimum is much deeper for 
this model. The fitted lines correspond to power laws of the form $ = a + bL w with exponents 
w = 0.97(3), w = 0.99(3), w = 0.99(3) and w = 0.92(3), respectively. The finest estimate to 
the expected critical exponent w = a/v = 1 is the one corresponding to the original definition 
(Eqs. 0-©) of the minimum subspaces. 
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FIG. 4: Relative deviations of WL ES schemes, with respect to the Metropolis algorithm defined 
in Eq. (|15|) . calculated from the order parameter at the critical temperature T c . The error bars used 
show the Metropolis relative uncertainties calculated as 5 standard deviations in the equilibrium 
regime. 
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FIG. 5: The same deviations, as in Fig. 0J again at T c , calculated now from the susceptibility. The 
error bars as in Fig. 01 show again the Metropolis relative uncertainties calculated as 5 standard 
deviations in the equilibrium regime. 
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FIG. 6: Susceptibility deviations at the corresponding pseudocritical temperatures. The error 
bars as in Figs. I4I5I 
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FIG. 7: Evolution of the estimates of the critical susceptibility with the WL iteration on a lattice 
of linear size L = 30. The filled (open) squares present the case mr-WLl (mr-WL2) which is 
the multi-range WL(l-24) (WL(12-24)) recipe, described in the text. Upper open triangles (mr- 
WL25) illustrate the evolution when a separation (<S) refinement of 16 spin-flips is applied between 
successive records in the (E,M) histograms. The last two cases, WL2S(tj) and WL25(3t,) (down 
filled triangles and open circles respectively), correspond to a simple one-range approach of different 
simulation times. The duration of each WL iteration was carefully chosen, so that saturation of 
the histogram fluctuation was well obeyed, as shown in the inset (j = 12). The solid line represents 
the estimates of xt c obtained by the Metropolis run, as explained in the text. 
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FIG. 8: 2D Ising model: Logarithmic scaling behavior of the finite-size extensions of the CRMES 
defined with the help of specific heat maxima according to the successive minimal approximations 
(Eq. (JSJ) and the analogous extensions of the CRMES defined with the help of the susceptibility 
maxima (Eq. 
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FIG. 9: 2D Ising model: Finite-size behavior of the extensions of critical minimum magnetic sub- 
spaces (CrMMS) obtained from the WL schemes, at the susceptibility pseudocritical temperatures, 
calculated with the help of the definition ()19|) and using r = 10~ 6 . The fitted line correspond to 
the power law 0.525 • L 175 . 
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FIG. 10: The same as in Fig. El at the exact critical temperature T c , including now for compar- 
ison the CrMMS corresponding to the order parameter probability distributions obtained by the 
Metropolis algorithm. The fitted line, used as a guide to the eye, is the power law 0.55 • L 175 , 
obtained by fixing the exponent to 1.75 and fitting the data L = 20 - 140 of the WL(N-fold: 12-24) 
scheme. Note the decline of the Metropolis CrMMS. 
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FIG. 11: Illustration of the very slow equilibration of the Metropolis algorithm in the tail regime. 
Time-development of the extensions of CrMMS (corresponding to r = 10~ 4 and r = 1CT 6 ) of the 
Metropolis algorithm. To observe the distance from the true- or tail-equilibrium we have divided 
the time-developing Metropolis extensions with the corresponding extensions of the WL(N-fold:12- 
24) scheme so that true-equilibrium occurs at the value 1. It is obvious that for small r and large 
L true-equilibrium of the Metropolis method is not attained, even for very long runs. 
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FIG. 12: The same as in Fig. ^2f° r L = 70 for a longer run. The dotted line shows the 
value obtained for r = 10~ 6 for the first 300 time-steps (compare with Fig. 1111) . 
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